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Abstract 

A compact finite-difference approximation to the unsteady Navier-Stokes 
equations in velocity-vorticity variables is used to numerically simulate a 
number of flows* These include two-dimensional laminar flow of a vortex 
evolving over a flat plate with an embedded cavity, the unsteady flow over an 
elliptic cylinder, and aspects of the transient dynamics of the flow over a 
rearward facing step. The methodology required to extend the two-dimensional 
formulation to three-dimensions is presented. 
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Introduction 


We present a solution method for the unsteady Navier-Stokes equations 
expressed in velocity-vorticity variables. Although this formulation is 
somewhat unusual, incompressible flow fields can be regarded as a realization 
of vorticity dynamics and as being driven by the production of vorticity at 
the boundaries. 

Most algorithms for the incompressible Navier-Stokes equations are couched 
in either stream function, vorticity variables in two-dimensions, or in so- 
called primitive variables of pressure and velocity in either two— or three- 
dimensions. When the stream function, vorticity formulation is used, accuracy 
can be degraded in differencing the stream function to obtain boundary values 
for the vorticity. Difficulties can arise, when using the primitive variable 
approach, in calculating the pressure boundary conditions and in maintaining 
the incompressibility condition. 

The formulation of the Navier-Stokes equations in terms of velocity and 
vorticity is an alternate approach. Dennis, Ingram, and Cook (1979) and Fasel 
(1980) have also used this formulation as the basis of numerical 
calculations. Dennis et al treated the steady-state problem in three 
dimensions and Fasel treated the time-dependent problem in two dimensions. In 
both of these studies, Poisson equations for the velocity components were 
derived from the kinematic definition of vorticity and used in the solution 
algorithm. In the numerical method developed by Gatski, Grosch, and Rose 
(1982) (GGR, hereafter) and used here, the kinematic definition of vorticity 
is used directly, along with the incompressibility condition of a divergence 
f ree-velocity field. These equations, coupled with the transport equation for 
the vorticity, form the basis of the algorithm. 
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The f inite-dif ference schemes used are members of a class of compact 
finite-difference schemes described by Keller (1974), Rose (1981), Philips and 
Rose (1982). They have been used by Wornom (1977), GGR, and Malik, Chuang, 
and Hussaini (1982). These compact schemes appear to have some advantages 
over more conventional central and upwind-difference schemes. Some of these 
advantages are that they are second-order accurate, independent of the value 
of the local cell Reynolds number, and require only values of the dependent 
variables in, and on, the boundaries of a single cell, thus easing the 
calculation of boundary conditions and the use of stretched grids. 

In addition to the formulation and test of the method (GGR), the algorithm 
has been applied to the laminar flow over an embedded cavity (Gatski and 
Grosch, 1984) and the unsteady forced flow of a shear layer (Mclnville, 
Gatski, and Hassan, 1984). Results presented here include the interaction of 
a Stuart vortex (Gatski, 1983) with an embedded cavity in laminar boundary 
layer flow, the unsteady flow past an elliptic cylinder, and the flow past a 
rearward facing step in a channel. 


Solution Method 

The basic development and formalism for the two-dimensional solution 
method are described in GGR. In the present work, the boundary specification 
has been altered in the vortex over an embedded cavity example to account for 
the incoming vortical motion, and in the elliptic cylinder example, the 
solution method is applied to the flow in a non-Cartesian coordinate system. 

For completeness, it is appropriate to give a brief outline of the method 
used in the solution of the two-dimensional problem. The governing 
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defferential equations for the solution set are those for an incompressible 
laminar flow, 


= o 


i ijk 

? = e Vj 


( 1 ) 


( 2 ) 


= Re 1 g jU ?^j k 


(3) 


where in Cartesian coordinates for example 


c 1 = (x,y,z); u 1 = (u,v,w); s 1 = (s 1 ,? 2 ,? 3 ). 


(4a,b,c) 


Re is a suitably defined Reynolds number based on the characteristic velocity 

and length scales of the flow and the kinematic viscosity of the fluid* 

Equations (1) through (3) are written in generalized tensor notation, with 

g^k the metric tensor, to expedite the coordinate transformation discussion 

to be presented later. In the following presentation of the numerical method, 

the discussion is limited to the two-dimensional Cartesian case and the x 

and y component velocities given by u and v, respectively, and the 
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nonzero vorticity component given by r, = £ . The finite-difference 
approximations to Equations (1) through (3) are given by GGR 

6 u 1 } + 6 v? = 0; 6 - 6 u? = (5a, b) 

x y x y 


with auxiliary conditions 
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(y -y )u? = 0; (y-y)v? = 0 (6a, b) 

x y x y 

and 

(6 t + ^x u "> 5 x + v ?)6 y )c " = Re *(6^” + 6 y i|>") (7a) 

6 = (y —1/2 Axq 6 H?; 6 c? = (y -1/2 Ayq 6 )i|>? (7b, c) 

x x xx y y y y 

with auxiliary conditions 


V 


n 


v- = 


v 


( 8 ) 


The dot subscript notations used in Equations (5) through (8) implies average 
and difference operations about cell centers, and the variables q x and q y 
are functions of the respective cell Reynolds numbers and are given by 


q = coth 0 - 0 ^ ; q = coth 0 - 0 ^ (9a, b) 

M x x x ’ y y y 

where 

0 = u.AxRe/2; 0 = v.AyRe/2 (10a, b) 

x y 

The solution method consists of an alternation between a velocity solver, 
described in GGR, for the velocity equations, Equations (5) and (6); and a 
vorticity solver, also described in GGR, for the vorticity equations, 

Equations (7) and (8). This methodology along with the appropriate time- 
dependent inflow boundary conditions described later constitute the technique 
used in solving the evolving vortex flow example. 

The solution method just described can also be directly applied to the 
solution of the unsteady flow over an elliptic cylinder by first transforming 
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the coordinate basis set from Cartesian to elliptic cylindrical. Consider the 
transformation 

x = (a^ - cosh £ cos n; y = (a^ - b^^^sinh £ sin n (lla,b) 

where a and b are the semi-major and semi-minor axes of the ellipse, 
respectively; and lines of constant 5 and n are confocal ellipses and 
hyperbolas, respectively. This system is particularly advantageous since the 
non-zero elements of the metric tensor are equal, that is 



( (a^ - b^)(sinh^ £ + sin^n) 
0 ° 


0 

(a^ - b^)(sinh^ £ 
0 


• 2 ^ 
sm n ) 


( 12 ) 


Thus the vorticity transport equation, expressed in $ > ti variables as well 
as the respective physical component velocities, is of the same form as the 
corresponding Cartesian vorticity equation. This then allows for the direct 
application of the difference equations, given in Equations (5) through (8) 
for the Cartesian system, to the differential equations governing the elliptic 
cylinder problem. 

Stream function and pressure values can be calculated from the computed 
velocity and vorticity variables. The stream function values are obtained 
directly from the velocity field by a simple integration of the velocity along 
the edge of each cell. This then determines the stream function values at the 
corners of each computational cell. The pressure field is more complex in 
that the pressure values are determined from the x and y momentum 
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equations and as such must be defined at the midpoint of the cell edges to be 
consistent with the velocity and voricity variable locations, A consistent 
means, therefore, of obtaining the pressure field is to use a solution 
procedure which is analogous to the one used in the velocity solver. This can 
be done by first writ ting the u and v momentum equations, or their non- 
Cartesian analogs in the case of the elliptic cylinder, as 

3p/3x = — (u + uu + vu + Re ) = L (u ) (13a) 
txy y x 

3p/3y = -(v + uv + vv - Re *£ ) = L (v,£) , (13b) 
txy x y 

where all the terras on the right side of Equations (13a) and (13b) are known 
at each time step. Next, recombine these equations into the more adaptable 
form 

3p/3x + 3p/3y = L x + ; 3p/3x - 3p/3y = (14a, b) 

Recall that the pressure variables are defined along the edges of each 
computational cell. If the pressure variables defined along the vertical 
edges (lines of constant x) are designated Py and those defined along the 
horizontal edges (lines of constant y) are designated P H» then the following 
discretization, of Equation (14) holds 


6 P + <$ P = J2(L + L ) ; 6 P - 6 P = ft(L 
xV yH x y xV yH x 


along with the auxiliary conditions 
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M x P V - Vh = 0 (15c) 

Since Equation (15) is an analog of Equations (5) and (6), the solver which 
was applied to the discretized velocity equations can be applied directly to 
the above discretized pressure equations. 

The results presented in the next section and from those presented in 
previous studies (Gatski and Grosch, 1984; Mclnville, Gatski, and Hassan, 
1984) show that the algorithm is adaptable to a rather wide class of two- 
dimensional unsteady flows. It is desirable to extend this methodology to 
three-dimensional unsteady flows. As was the case in the two-dimensional 
problem, the differential equation set, Equations (1) through (3), are used 
directly in the three-dimensional problem. Equations (1) and (2), which once 
again constitute the velocity solver, are discretized using box— variables to 
represent the velocities and assigning the vorticity values to the center of 
each computational box. Such a formulation produces an over— determined 
system; however, Fix and Rose (1984) have shown that such a finite-difference 
approximation yields a least squares solution which is second order 
accurate. This velocity field is then used to calculate vorticity boundary 
conditions through second order accurate difference approximations at the 
boundaries. Before the vorticity transport equations can be brought into the 
solution sequence a modification to the form of the equations needs to be 
made. This is necessary because the three-dimensional vorticity transport 
equations contain a vortex stretching term which negates the direct use of the 
finite-difference basis set which was used in the two-dimensional transport 
equation. Recall, however, that as the vorticity equations are solved over a 
time step At, the form of the equations allows for the introduction of an 
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integrating factor (Rose, private communication) which transforms the 
vorticity transport equations into simple advection-diffusion equations. This 
system can then be solved for the three component vorticities in a completely 
analogous manner to the simple two-dimensional equation. 

This completes the outline of the methodology that can be used in the 
solution of a class of two- and three-dimensional unsteady flows. In the next 
section, the application of this solution technique to a few relevant flow 
problems is presented. 


Computational Results 

In this section some results from the numerical solution of three two- 
dimensional unsteady flows will be presented. The two-dimensional flows 
include the evolution of a vortical structure over an embedded cavity, the 
unsteady flow over an elliptic cylinder, and some aspects of the flow over a 
rearward facing step in a channel. 

Consider first the evolution of a vortical structure over an embedded 
cavity. Such a flow is of interest since it serves as a qualitative model of 
a large scale vortical structure, omnipresent in wall bounded turbulent shear 
flows, evolving over an isolated surface roughness represented by an embedded 
cavity. The reference flow field, in the absence of the perturbing vortical 
structure, is a laminar boundary layer flow over an embedded cavity. This 
flow has been extensively studied numerically in Gatski and Grosch (1984) 
where such flow field characteristics as wall shear stress and pressure 
distributions as well as the usual vorticity and stream function contours have 
been presented. In the present study, a Stuart vortex is introduced at the 
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inflow boundary in a consistent mathematical manner through the method of 
matched asymptotic expansions (Gatski, 1983); and flow characteristics similar 
to those obtained in the unperturbed case are obtained at different times in 
the evolution of the vortex over the cavity. Figure 1 shows the vorticity, 
stream function and pressure contours of the flow when the vortex has evolved 
to a point directly above the embedded cavity. Here, the cavity is square and 
has a depth of one-half of the inflow boundary layer thickness. The vorticity 
contours are shown in Figure la. The figure shows the main vortical motion 
above the cavity as well as a remnant of an induced vortical region downstream 
of the main motion. This induced structure was originally formed by the 
interaction of the Stuart vortex with the bounding wall; however, as the main 
vortex evolved to the point shown in the figure, the strength of this induced 
vortex was weakened by the embedded cavity. Figure lb shows an enlarged view 
of the motion in the cavity, as represented by the stream function contours, 
when the vortex is in the position shown in Figure la. It is, of course, 
realized that any interpretation of the unsteady motion of the fluid in terms 
of the stream function must be done with caution; however, to be consistent 
with the previous work done on this same flow geometry (Gatski and Grosch, 
1984) but without the vortex, the stream function contours are qualitatively 
informative. In the unperturbed flow, the cavity vortex was bounded by the 
cavity walls and by a slightly convex zero streamline, which indicated a local 
acceleration of the flow in the vicinity of the cavity. The results shown in 
Figure lb indicate that the vortex in the boundary layer causes the cavity 
vortex to lift out of the cavity. This type of vortex action is a significant 
factor in the alteration of the drag characteristics of turbulent flow over 
surface roughnesses. In Figure lc is shown the pressure contours in the 
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boundary layer. Ihe figure clearly shows the low pressure regions associated 
with the area beneath the vortices and the high pressure region associated 
with the flow which is directed toward the wall by the main vortex. This 
pictoral representation of the perturbed cavity flow is intended to serve as 
an indicator of the rather complex dynamics which takes place in this type of 
flow field. The additional computational results which have been obtained 
from this study have allowed for the accurate calculation of such quantities 
as wall shear stress and wall pressure distributions which are used to 
quantify the drag characteristics of such flows. 

A second example is the impulsive start of the flow over an elliptic 
cylinder. This is an example of the unsteady separation of an external 
flow. The impulsive start of a circular cylinder is, of course, a classic 
problem, but there seems to be no numerical work for the elliptic cylinder. 
We are beginning a systematic program of such calculations, in the course of 
which we intend to vary the Reynolds number, slenderness of the ellipse, and 
angle of attack. 

Some preliminary results are shown in Figure 2 at a time shortly after the 
beginning of separation. Here the Reynolds number, with the length scale 
based on the semi-major axis, is 100, the ratio of major to minor axis is 2, 
and the angle of attack is zero. The computational domain is 0 < n < 2tt , 
0.5493 < E, < 4.2493; the outer boundary is approximately that of a circle with 
radius = 22a. There are 120 cells in the n direction with a width of 
ir/60, and 50 cells in the % direction. Variable cell heights, ranging from 
0.01 to 1.0, are used in the £ direction in order that the boundary layer 
can be resolved. The velocity and vorticity are specified on the forward 
portion (it/ 2 < q < 3it/2) of the outer boundary. Outflow boundary conditions 
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are applied to the remainder of the outer boundary. The initial condition is 
that the flow is that given by the potential solution. 

The vorticity distribution, as shown in Figure 2a, is determined by a 
balance between diffusion and advection. Vorticity is produced at the 
boundary and is diffused away. Simultaneously it is advected towards the rear 
of the ellipse. This leads to the appearance of a rather thin region of non- 
zero vorticity on the forward part of the ellipse and the long "tails" of 
vorticity streaming from the body. This has been observed in both start-up 
and steady flow past circular cylinders and is discussed by Batchelor 
(1967). In contrast, the stream function results given in Figure 2b indicate 
that the overall velocity field is only slightly different from that of the 
potential flow. A thin region of separation is also apparent in Figure 2, but 
at this early time a viscous wake has not yet formed. Finally, it can be seen 
from the results shown in this figure that the flow field is nearly, but not 
quite symmetric about ri =0, tt . Symmetry conditions are not imposed on the 
solution, but the velocity and vorticity are required to be periodic in n 
and they are symmetric at t = 0. This asymmetry is due to small 
perturbations which are caused by the sweep direction bias in the solvers. 

Another example of practical relevance is the flow in a channel with a 
backward facing step. This is a simple prototype of a separating internal 
flow. Figures 3a and 3b show the steady-state distribution of the stream 
function and vorticity contours for this flow at Re = 300. The Reynolds 
number, Re, length scale is the height of the inflow channel and the velocity 
scale is the maximum speed at inflow. The inflow boundary conditions are a 
parabolic velocity profile and a linear distribution for the vorticity. The 
outflow boundary condition assumed a zero cross-stream component of velocity, 
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v, and a simple advective flux of vorticity across the outflow. The geometry 
of Figure 3 has been considerably distorted in order to display the results 
clearly. The height of the inflow channel and of the step is one unit. The 
length of the inflow region, upstream of the step, is four units, and the 
length of the channel downstream of the step is 30 units. The computational 
cells are of constant height, 0.04, in the cross-channel direction so that 
there are, across the channel, 25 cells in the entrance region and 50 cells 
downstream of the step. The importance of having good resolution near the 
step and the reattachment point, combined with the great length (34 units) of 
the computational domain required that variable grids be used along the 
channel. A total of 130 cells were used with the cell width varying from 0.01 
to 0.5. 

As can be seen from the results of Figure 3, there is, at this Reynolds 
number, virtually no upstream effect of the step. It was found that the 
deviation from a parabolic velocity and linear vorticity profile was less than 
10 -3 one unit upstream of the step. Thus, an entrance region four units long 
is more than adequate. In contrast, a computational domain which is 30 units 
long downstream of the step is, at this Re, only just adequate. Originally, 
the region downstream of the step was taken to be 20 units long. It was found 
that this was too short. The velocity and vorticity profiles in the outflow 
region were distorted and not symmetric about the centerline of the channel. 
Considerable numerical experimentation led to the conclusion that this was due 
to the channel being too short; thus the computational domain was lengthened 
by 10 units. 

The present results show that the flow over the step generates a rather 
weak corner vortex, at least at this low Re. The speed of the reverse flow 
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in this vortex is less than 10 percent of the maximum inflow speed. It 
appears that there is a diffusive dominance in the cross-channel direction; and 
an advective dominance along the channel. One result of this is that the 
effect of the step persists, in the lower half of the channel for very large 
distances downstream. In the upper half of the channel, there is a region of 
decelerating flow downstream of the step and over the recirculation zone. 
This deceleration, combined with the frictional drag of the upper wall gives a 
region of near separation on the upper wall. We have found that at slightly 
higher Re, a zone of laminar separation and reattachment forms on the upper 
wall. Beyond the separation or near separation zone, the readjustment of the 
flow is due to a cross-channel diffusion of strearawise momentum. A long 
length of channel is required for this process to be completed. 

As has been shown, this flow field is stable, but that does not prevent 
the formation of stable, that is decaying, shear waves. Any impulsive change 
in the flow will excite these transient waves. An example is shown in Figure 
4. The flow field at Re = 300 was perturbed by decreasing the viscosity so 
that the Reynolds number was instantaneously changed from 300 to 500. 
Contours of the instantaneous values of stream function and vorticity are 
given in Figures 4a and 4b, respectively. These results are at approximately 
20 time units after the viscosity was changed and show large amplitude shear 
waves in the channel. Note that these waves are downstream of the separation 
zone behind the step; which is, in fact, where they were formed. A comparison 
of Figures 4a and 4b show that the vorticity field is a more sensitive 
indicator than the stream function. The waves near the downstream boundary 
are quite visible in the vorticity contours, while they are nearly absent in 
the same region of the plot of the instantaneous streamlines. There is a 
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maximum of approximately 0.73 in the stream function downstream of the step 
and near the upper boundary, which is considerably above 2/3, the maximum 
value of the stream function at inflow and outflow. Finally, there is a 
separation region on the upper wall, as shown by the region of negative 
vorticity, which is assocsiated with this stream function maximum. 

Finally, the three-dimensional formulation which was described in the 
previous section is presently being applied to the study of the time-dependent 
behavior of a single Taylor-Green vortex. As is known such a vortex is 
described by a simple diffusion equation for each vorticity component and 
serves as an excellent check on the diffusive characteristics of the numerical 
algorithm in three dimensions. In addition, since the structure of the 
Taylor-Green vortex is known analytically, this model problem is also a useful 
means of checking the accuracy of the overall solution method. Work on the 
Taylor-Green problem is in its early stages and will be presented elsewhere. 


Concluding Remarks 

The results presented in this study have shown that the compact vorticity- 
velocity difference formulation used is readily adaptable to a rather large 
class of unsteady fluid flow problems. In addition, the results from this and 
previous studies have shown that the computed results are of sufficient 
accuracy to describe detailed dynamic features of these flows. Future efforts 
in this area include further development of the three-dimensional formulation 
and methodology, and the adaptation of this compact difference formulation to 
the new generation of concurrent processors with the goal of extending the 
range of real flows that can be treated in the present context. 
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Figure 2: (2a) Vorticity contours in near 

field of elliptic cylinder (contour levels 
-6.40 to 6.40); (2b) Stream function con- 
tours in near field of elliptic cylinder 
(contour levels -0.56 to 0.56). 


Figure 1: (la) Vorticity contours in boundary 

layer (contour levels -1.8 to 0.0); (lb) Stream 
function contours in embedded cavity (contour 
levels -0.002 to 0.005); (lc) Pressure contours 
in boundary layer (contour levels -0.017 to 
M0). 
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